#!/usr/bin/env python # coding: utf-8 # In[1]: import copy maxima_calculus('algebraic: true;') # This helps Sage to deal with fractions phi = (sqrt(5)+1)/2 # In[2]: # Quaternion multiplication def quat_mult(z, w): (a,b,c,d,e,f,g,h) = (z[0], z[1], z[2], z[3], w[0], w[1], w[2], w[3]) A = a*e - b*f - c*g - d*h B = a*f + b*e + c*h - d*g C = a*g + c*e + d*f - b*h D = a*h + d*e + b*g - c*f return [A,B,C,D] # Simplify expressions for quaternions (to be used only below - it does not work in general) def simplify_quaternion(q): p = q.copy() for i in range(4): b = q[i] if b != 0 and b != 1/2 and b!= -1/2: c = b.full_simplify() p[i] = c return p # In[3]: # Construct the set S of 6 quaternions in the binary icosahedral group close to 1 S = [[], [], [], [], [], []] S[0] = [phi/2, 0, 1/(2*phi), 1/2] S[1] = [phi/2, 0, 1/(2*phi), -1/2] S[2] = [phi/2, 1/(2*phi), 1/2, 0] S[3] = [phi/2, 1/(2*phi), -1/2, 0] S[4] = [phi/2, 1/2, 0, 1/(2*phi)] S[5] = [phi/2, -1/2, 0, 1/(2*phi)] for i in range(6): S[i] = simplify_quaternion(S[i]) for q in S: print (q) # Construct the binary icosahedral group of order 120, that is generated by S I = [] for q in S: I.append(q.copy()) sono_tutti = False while sono_tutti == False: for p in I: for q in I: r = quat_mult(p,q) r = simplify_quaternion(r) if I.count(r) == 0: I.append(r) length = len(I) print (length, r) if length == 120: sono_tutti = True break if sono_tutti == True: break # Sort the facets of the 120-cells following levels I.sort() print(I) # In[4]: # Determine the adjacency matrix of the 120-cell I # G[i][j] = 0 or 1 depending on whether the i-th and j-th vertices of I are close G = matrix(120,120) for i in range(120): print (i, end = ' ') for j in range(i+1,120): product = sum([I[i][t]*I[j][t] for t in range(4)]) if product == 1/4*sqrt(5) + 1/4: G[i,j] = 1 G[j,i] = 1 # In[4]: # Determine the 14.400 isometries of the 120-cell isometries = get_isometries(G) # In[9]: # Store the isometries on a file import pickle f = open("Isometries", "wb") pickle.dump(isometries, f) # In[58]: # Load the isometries from the file import pickle g = open("Isometries", "rb") # isometries_copy = pickle.load(g) isometries = pickle.load(g) print(len(isometries)) # In[59]: # Construct all the faces of the 120-cell, as the dual of the flag complex of G # faces[i] = (4-i)-uples of numbers in [0, ..., 119] faces = [[], [], [], []] faces[3] = [i for i in range(120)] for i in range(120): for j in range(i+1,120): if G[i,j] == 1: faces[2].append([i,j]) for i in range(120): for j in range(i+1,120): for k in range(j+1,120): if G[i,j] == 1 and G[i,k] == 1 and G[j,k] == 1: faces[1].append([i,j,k]) for i in range(120): for j in range(i+1,120): for k in range(j+1,120): for l in range(k+1,120): if G[i,j] == 1 and G[i,k] == 1 and G[j,k] == 1 and G[i,l] == 1 and G[j,l] == 1 and G[k,l] == 1: faces[0].append([i,j,k,l]) print (len(faces[0]), len(faces[1]), len(faces[2]), len(faces[3])) # In[11]: # Identify the orbits of the faces after the identification in the Davis manifold # Warning: uses the fact that facets i and 119-i are opposite # The algorithm is quite slow orbits = [[], [], [], []] # Orbits of facets for i in range(120): for j in range(i+1, 120): if I[i] == quat_mult(I[j], [-1,0,0,0]): orbits[3].append([i,j]) # Orbits of pentagons centers = [] for f in faces[2]: new_center = [(I[f[0]][i] + I[f[1]][i])/2 for i in range(4)] centers.append(new_center) for i in range(len(faces[2])): center = centers[i] adjacent_facets = faces[2][i] gia_trovato = False for o in orbits[2]: for f in o: if set(f) == set(adjacent_facets): gia_trovato = True if gia_trovato == False: new_orbit_centers = [center[:]] keep_going = True print (i, center, adjacent_facets) while keep_going: keep_going = False for c in new_orbit_centers: for t in range(2): base_vector = I[adjacent_facets[t]] Fourier = sum([c[i] * base_vector[i] for i in range(4)]) # base_vector is unitary new_c = [c[i] - 2 * Fourier * base_vector[i] for i in range(4)] new_c = simplify_quaternion(new_c) if new_orbit_centers.count(new_c) == 0: new_orbit_centers.append(new_c) keep_going = True print (new_c) new_orbit = [] for j in range(len(faces[2])): if new_orbit_centers.count(centers[j]) == 1: new_orbit.append(faces[2][j]) print("Orbita:", new_orbit) orbits[2].append(deepcopy(new_orbit)) # Orbits of edges centers = [] for f in faces[1]: new_center = [(I[f[0]][i] + I[f[1]][i] + I[f[2]][i])/2 for i in range(4)] centers.append(new_center) for i in range(len(faces[1])): center = centers[i] adjacent_facets = faces[1][i] gia_trovato = False for o in orbits[1]: for f in o: if set(f) == set(adjacent_facets): gia_trovato = True if gia_trovato == False: new_orbit_centers = [center[:]] keep_going = True print (i, center, adjacent_facets) while keep_going: keep_going = False for c in new_orbit_centers: for t in range(3): base_vector = I[adjacent_facets[t]] Fourier = sum([c[i] * base_vector[i] for i in range(4)]) # base_vector is unitary new_c = [c[i] - 2 * Fourier * base_vector[i] for i in range(4)] new_c = simplify_quaternion(new_c) if new_orbit_centers.count(new_c) == 0: new_orbit_centers.append(new_c) keep_going = True print (new_c) new_orbit = [] for j in range(len(faces[1])): if new_orbit_centers.count(centers[j]) == 1: new_orbit.append(faces[1][j]) print("Orbita:", new_orbit) orbits[1].append(deepcopy(new_orbit)) # In[13]: # We get 60 facets, 144 pentagons, and 60 edges. print("Facets:") print(orbits[3]) print("Pentagons:") print(orbits[2]) print("Edges:") print(orbits[1]) # In[14]: # Store the cellularization on a file import pickle f = open("Orbits", "wb") pickle.dump(orbits, f) # In[60]: # Load the cellularization from the file import pickle g = open("Orbits", "rb") # orbits_copy = pickle.load(g) orbits = pickle.load(g) # In[61]: # For each of the 144 pentagons, write the 5 dodecahedra that contain it # and the 5 edges that are contained in it. # Note that for each pentagon there is the opposite one with the same data poset = [] for i in range(144): o = orbits[2][i] p = [[], []] for j in range(60): # Determina i dodecaedri che contengono il pentagono trovato = False for t in range(2): for u in range(5): if orbits[2][i][u].count(orbits[3][j][t]) == 1: p[0].append(j) trovato = True break if trovato == True: break for j in range(60): # Determina gli spigoli contenuti nel pentagono trovato = False for t in range(20): contenuto = True for l in range(2): if orbits[1][j][t].count(orbits[2][i][0][l]) == 0: contenuto = False if contenuto == True: p[1].append(j) trovato = True break poset.append(deepcopy(p)) print (poset) # In[9]: # Recursive function that construct all the 14400 inside-out isometries, as permutations p of [0, ..., 59] # The dodecahedron i is sent to the edge p[i] def inside_out(level = 0, isometry = [], isometries = [], n = 0): if n == 0: n = 60 isometry = [0 for i in range(n)] isometries = [] if level == n: isometries.append(isometry[:]) print (len(isometries), isometry) for i in range(n): if isometry[:level].count(i) == 0: isometry[level] = i print (isometry, end='\r', flush=True) # ACHTUNG is_ok = True for f in range(144): parents = [] for j in range(5): if poset[f][0][j] <= level: parents.append(poset[f][0][j]) trovata = False for g in range(144): contains = True for parent in parents: if poset[g][1].count(isometry[parent]) == 0: contains = False if contains == True: trovata = True break if trovata == False: is_ok = False break if is_ok == True: inside_out(level + 1, isometry, isometries, n) if level == 0: return isometries # In[10]: # Call the recursive function inside_out() # In[23]: # I write here the first inside-out isometry the routine finds above # So there is no need to run it again IO = [0, 1, 5, 15, 46, 53, 54, 52, 51, 45, 16, 6, 2, 55, 50, 44, 43, 19, 25, 9, 10, 26, 20, 23, 29, 13, 12, 28, 22, 39, 42, 49, 56, 3, 7, 17, 32, 31, 33, 36, 34, 35, 18, 8, 4, 14, 24, 30, 40, 37, 47, 48, 57, 58, 38, 41, 21, 27, 59, 11] print(IO) # In[28]: # Construct 72 oriented great circles as oriented cosets # of the 6 initial order-10 cyclic subgroups # generated by the elements in S. # Each subgroup has 12 cosets # C[12*i+j] is the j-th coset of the i-th subgroup, with i = 0,...,5 and j = 0,...,11 # # This list could also be deduced more easily from orbits[2], which was written later. # We keep it as it is just in case some parts could be useful. C = [] for i in range(6): # Construct the order-10 cyclic subgroup H generated by q = S[i] q = S[i] p = [1,0,0,0] H = [] for i in range(10): p = quat_mult(p,q) p = simplify_quaternion(p) H.append(p) # Construct the 12 cosets of H CH = [] for p in I: pH = [] for q in H: pq = quat_mult(p,q) pq = simplify_quaternion(pq) pH.append(pq) already_there = False for coset in CH: for a in pH: for b in coset: if a == b: already_there = True if already_there == False: CH.append(pH.copy()) C.append(pH.copy()) print (len(C), pH) if len(CH) == 12: break # Sort the 72 circles individually so that each circle starts from a minimmum dodecahedron # Then sort the 72 circles for c in C: minimo = 0 for i in range(10): if c[i] < c[minimo]: minimo = i d = [c[(minimo + t) % 10] for t in range(10)] for t in range(10): c[t] = d[t] C.sort() # Translate C so that the elements of I are numbers from 0 a 119 # The new list is Ccopy Ccopy = [] for c in C: cnew = [I.index(c[j]) for j in range(10)] Ccopy.append(deepcopy(cnew)) print (cnew) # In[16]: # Store Ccopy on a file import pickle f = open("Ccopy", "wb") pickle.dump(Ccopy, f) # In[62]: # Load Ccopy from the file import pickle g = open("Ccopy", "rb") # Ccopycopy = pickle.load(g) Ccopy = pickle.load(g) # In[23]: # Determine the intersection form Q in the basis of the great circles. # The intersection number of two circles is 1, 0, or -1: # it is 0 if they intersect, # It is 1 or -1 if they are disjoint, and in this case it is their algebraic intersection number. M = [([0]*72) for i in range(72)] # M is a 72x72 matrix of zeroes for i in range(72): for j in range(72): a, b = C[i], C[j] A = matrix([a[0], a[1], b[0], b[1]]) det = A.determinant() M[i][j] = det if M[i][j] != 0: if M[i][j] > 0: M[i][j] = 1 if M[i][j] < 0: M[i][j] = -1 print (i, end = ' ') # We store this information in an integer matrix Q Q = Matrix(ZZ, M) rank(Q) # We check that det(Q) = 2^32 d = Q.determinant() print (d, d == 2^72) # In[24]: # Store the intersection form Q on a file import pickle f = open("Q", "wb") pickle.dump(Q, f) # In[63]: # Load the intersection form Q from the file import pickle g = open("Q", "rb") # Q_copy = pickle.load(g) Q = pickle.load(g) # In[68]: # Check that the determinant is 2^72 Q.determinant(), 2^72 # In[7]: # Determine the isomorphisms of C induced by the 14.400 isometries of the 120-cell # and by the reversal of signs (the matrix -I) # Cisom[i] = [[7,-1], [98,1], [10,-1], ... ] # indicates that the i-th isomorphism sends the circle 0 to the circle 7 inverting its orientation, etc. Cisom = [] counter = 0 for p in isometries: new_isom = [] for i in range(72): c = Ccopy[i] cimage = [p[c[j]] for j in range(10)] for j in range(72): they_are_equal = True for k in range(10): if Ccopy[j].count(cimage[k]) == 0: they_are_equal = False break if they_are_equal == True: cambio_orientazione = 1 for k in range(10): if cimage[k] == Ccopy[j][0] and cimage[(k+1) % 10] == Ccopy[j][9]: cambio_orientazione = -1 new_isom.append([j,cambio_orientazione]) break print (counter, end='\r', flush=True) counter = counter + 1 Cisom.append(deepcopy(new_isom)) # We sort Cisom. # There are duplicates. Change the duplicates by inverting signs: this is allowed since all inequalities # form a symmetric set under all signs inversion. So we get 14.400 isomorphisms overall. Cisom.sort() for i in range(14399,0,-2): for j in range(72): Cisom[i][j][1] = - Cisom[i][j][1] print(len(Cisom)) # Add to each isomorphism an initial flag # that indicates the list of the i such that if p[k]=i then p[j] i: is_good = False break if p[j][0] == i: break if is_good: flags.append(i) new_p = [flags[:], p] CCisom.append(deepcopy(new_p)) # print(new_p) # Finally, construct the list Isom_efficient, which will be used by the algorithm # Isom_efficient[i] contains the list of isomorphisms whose initial flag contains i # # We avoid constructing Isom_efficient[71] because it would contain all the 144000 isometries # It takes a lot of time to construct it and it would be useless anyway # Hence Isom_efficient[71] is an empty list Isom_efficient = [] for i in range(72): new_list = [] if i != 71: # This condition ensures that Isom_efficient[71] is empty (see above) for p in CCisom: if p[0].count(i) == 1: new_list.append(deepcopy(p[1])) Isom_efficient.append(deepcopy(new_list)) print (i, end = ' ') # Removes the identity that is useless for our purposes for i in range(71): Isom_efficient[i].remove([[0, 1], [1, 1], [2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1], [8, 1], [9, 1], [10, 1], [11, 1], [12, 1], [13, 1], [14, 1], [15, 1], [16, 1], [17, 1], [18, 1], [19, 1], [20, 1], [21, 1], [22, 1], [23, 1], [24, 1], [25, 1], [26, 1], [27, 1], [28, 1], [29, 1], [30, 1], [31, 1], [32, 1], [33, 1], [34, 1], [35, 1], [36, 1], [37, 1], [38, 1], [39, 1], [40, 1], [41, 1], [42, 1], [43, 1], [44, 1], [45, 1], [46, 1], [47, 1], [48, 1], [49, 1], [50, 1], [51, 1], [52, 1], [53, 1], [54, 1], [55, 1], [56, 1], [57, 1], [58, 1], [59, 1], [60, 1], [61, 1], [62, 1], [63, 1], [64, 1], [65, 1], [66, 1], [67, 1], [68, 1], [69, 1], [70, 1], [71, 1]]) # In[15]: # Truncate the list Isom_efficient to a cutoff, otherwise it is too long to be written on a file Cutoff = 40 tot = 0 for i in range(Cutoff): tot = tot + len(Isom_efficient[i]) print (tot) for i in range(Cutoff, 72): Isom_efficient[i] = [] # In[19]: # Store Isom_efficient on a file import pickle f = open("Isom_efficient", "wb") pickle.dump(Isom_efficient, f) # In[20]: # Load Isom_efficient from the file import pickle g = open("Isom_efficient", "rb") # Isom_efficient_copy = pickle.load(g) Isom_efficient = pickle.load(g) # In[26]: # We determine the action R of the inside out isometry on homology, which is a 72 x 72 matrix # that represents it with respect to the basis of the oriented circles # Let OI be the inverse of IO OI = [0 for i in range(60)] for t in range(60): OI[IO[t]] = t print("OI =", OI) # The isometry IO acts on the 72 circles, because it sends a circle of 5 dodecahedra, determined by two oppposite # pentagons, into a set of 5 edges, that are also determined by two opposite pentagons, that in turn # determine a great circle. # We determine this isomorphism and call it "action" poset_semplificato = [] for i in range(144): new_poset = [[], []] for j in range(5): new_poset[0].append(orbits[3][poset[i][0][j]][0]) new_poset[0].append(orbits[3][poset[i][0][j]][1]) for j in range(5): new_poset[1].append(orbits[3][OI[poset[i][1][j]]][0]) new_poset[1].append(orbits[3][OI[poset[i][1][j]]][1]) poset_semplificato.append(deepcopy(new_poset)) # for p in poset_semplificato: # print (p) action = [0 for i in range(72)] for i in range(72): c = Ccopy[i] for p in poset_semplificato: if set(p[0]) == set(c): for j in range(72): if set(p[1]) == set(Ccopy[j]): action[i] = j print("Action =", action) # Determine the action R of the inside-out isometry on H2 Qinv = 2 * Q.inverse() # print (Qinv) R = Matrix(QQ, 72) for i in range(72): for j in range(72): R[i,j] = Qinv[i,action[j]] # The matrix R is not yet correct because we have not cared about signs. # So we need to change the signs of the colomns of R, and we do it now. for k in range(72): for i in range(72): x = R.column(k) * Q * R.column(i) print(x, Q[k,i]) if x != Q[k,i]: for j in range(72): R[j,i] = - R[j,i] # Finally we check that R is a Q-isometry: this should return "True" print(R.transpose() * Q * R == Q) # In[29]: # We determine the 360 inequalities produced by the elements b_i # Each of the 60 dodecahedra contains 6 large hexagonal paths that cut the dodecahedron in two equal pieces. # Each large hexagonal path represents a b_i, so they are 60x6 = 360 # Each b_i intersects only the 6 elements a_i that go through that dodecahedron, with a sign. # The matrix B is 360 x 72 and B[i][j] is the intersection between b_i and a_j B = [] # Beven = [] # Tiene conto solo della parità, viene una matrice 60 x 72 for g in I: if g > [0,0,0,0]: # Ne selezioniamo solo metà perché facce opposte sono identificate circles = [] # This list will contain the indices of the 6 circles that contain g for i in range(72): if C[i].count(g) == 1: circles.append(i) for j in range(6): # Each of the 6 circles determines a b_k c = C[circles[j]] # This is the circle that determines b_k cseq = c[(c.index(g) + 1) % 10] v = [0 for i in range(72)] # v[i] = for i in range(6): d = C[circles[i]] dseq = d[(d.index(g) + 1) % 10] c_delta = [cseq[i] - g[i] for i in range(4)] d_delta = [dseq[i] - g[i] for i in range(4)] product = sum([c_delta[i] * d_delta[i] for i in range(4)]) print (cseq, dseq, product) if product > 0: v[circles[i]] = 1 if product < 0: v[circles[i]] = -1 B.append(v.copy()) # if j == 0: # w = [abs(a) for a in v] # Beven.append(w.copy()) # We now construct BB and BBC from B. These will be useful to speed up the routine # BB[i] = [a, B[i]] # where a is the maximum index such that B[i][a] != 0 # BBC[a][j][i] = B[i] # where B is the j-th B whose maximum index is a BB = [] for b in B: a = 0 for i in range(72): if b[i] != 0: a = i new_b = [a, b[:]] BB.append(new_b) BB.sort() BBC = [[] for i in range(72)] for b in BB: new_vector = [] for j in range(72): if b[1][j] != 0: new_vector.append([j, b[1][j]]) # BBC[b[0]].append(deepcopy(b[1])) BBC[b[0]].append(deepcopy(new_vector)) # Find the inequalities that are dual to those of B via the inside-out map BBdual = [] for i in range(len(BB)): v = vector(BB[i][1]) w = R.transpose() * v print(w) BBdual.append(w) # In[30]: # Store BBC and BBdual on a file import pickle f = open("BBC", "wb") pickle.dump(BBC, f) f = open("BBdual", "wb") pickle.dump(BBdual, f) # In[31]: # Load BBC and BBDual from the file import pickle g = open("BBC", "rb") BBC_copy = pickle.load(g) # BBC = pickle.load(g) g = open("BBdual", "rb") BBdual_copy = pickle.load(g) # BBdual = pickle.load(g) # In[35]: # Recursive function that searches for even integral classes x = x_ia_i such that: # || <= 2, that is |x_i| <= 1 # || <= 2, which is calculated using B # >= M # # x even implies that x_i = -1,0,1 and is even, and we check that. # M = 0 # We set M = 0, but we could put some other value here def search_classes(level = 0, vec = [], n = 0, answer = [], X = [], mancanti = []): if n == 0: n = 72 answer = [] vec = [-10 for i in range(72)] X = [[[] for i in range(72)] for j in range(72+360)] # ACHTUNG mancanti = [[[] for i in range(72)] for j in range(72+360)] for i in range(72): for j in range(72): temp = [abs(Q[i][k]) for k in range(j+1,72)] mancanti[i][j] = temp.count(1) for i in range(360): for j in range(72): temp = [abs(BBdual[i][k]) for k in range(j+1,72)] mancanti[72+i][j] = sum(temp) if level < n: for t in range(-1,2,1): # for t in range(1,-2,-1): # To be used if you want to run the counter in the other direction print (vec, end='\r', flush=True) vec[level] = t for i in range(72): if level == 0: X[i][0] = Q[i][0] * vec[0] else: X[i][level] = X[i][level-1] + Q[i][level] * vec[level] for i in range(360): if level == 0: X[72+i][0] = BBdual[i][0] * vec[0] else: X[72+i][level] = X[72+i][level-1] + BBdual[i][level] * vec[level] is_good = True # Here we put the filters : # if level == 20: # To be used if you want the counter to start from a 'start' # start = [-1, -1, -1, -1, -1, -1, -1, -1, 0, -1, 0, -1, 1, 0, -1, 0, 0, -1, 1, 0, 1, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, -1, -1, 0, 0, 1, -1, 1, -1, 0, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -10] # if vec < start: # is_good = False for b in BBC[level]: product = sum([vec[b[i][0]] * b[i][1] for i in range(6)]) if product != -2 and product != 0 and product != 2: is_good = False break if is_good == True: for p in Isom_efficient[level]: j = 0 while j < 72 and p[j][0] <= level: a = vec[p[j][0]]*p[j][1] if a < vec[j]: is_good = False break if a > vec[j]: break j = j + 1 if is_good == False: break if is_good == True: for i in range(72+360): if abs(X[i][level]) - mancanti[i][level] > 2: is_good = False break if is_good == True: search_classes(level + 1, vec, n, answer, X, mancanti) if level == 0: return answer if level == n: is_good = True # Here put the final filters : w = vector(vec) for i in range(72): product = w*Q[i] if abs(product) > 2: is_good = False return for i in range(360): product = w*BBdual[i] if abs(product) > 2: is_good = False return product = w*Q*w if abs(product) < M: is_good = False return if is_good == True: answer.append(vec[:]) print ("found:", abs(product), vec) # In[37]: # This is the main routine L = search_classes()